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In  this  Report,  we  summarize  results  obtained  with  support  of  ONR  Contract  N00014-C- 
90-0039.  Two  general  classes  of  problems  are  studied.  First,  we  have  performed  a  variety 
of  numerical  simulations  of  turbulence  in  the  presence  of  a  free  surface.  Second,  we  have 
performed  Very  Large  Eddy  Simulations  of  the  turbulent  flow  past  prototype  propeller 
conflgurations. 

The  structure  of  turbulence  interacting  with  a  free  surface  is  an  important  problem  for 
a  variety  of  Navy  applications.  First,  the  free  surface  may  have  significant  effects  on  the 
structure  of  turbulence  and  thereby  effect  transport  properties  in  ocean  environment.  Such 
knowledge  is  important  to  understand  such  varied  characteristics  as  sound  propagation, 
transport  of  nutrients,  contaminants,  and  other  tracers  in  the  upper  ocean.  Second,  the 
turbulence  may  itself  affect  free  surface  characteristics.  Such  effects  may  influence  our 
ability  to  discern  sub-surface  phenomena  from  observations  of  the  surface  waves.  Such 
information  is  critical  for  the  proper  interpretation  of  sea-surface  observations  obtained  via 
various  remote  sensing  techniques.  Third,  the  basic  structure  of  free  surface  turbulence 
has  remained  a  hitherto  relatively  unexplored  area  in  which  advances  may  influence  our 
fundamental  understanding  of  turbulence  processes. 

In  Section  2  we  summarize  our  results  on  free  surface  turbulence  through  a  variety  of 
simulations  of  turbulent  channel  flows  with  a  free  surface. 
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In  Section  3,  we  survey  our  numerical  methods  for  solving  moving  boundary  problems 
using  spectral  elements  and  finite  volume  techniques.  In  particular,  we  describe  the  RNG 
K-e  model  for  treatment  of  transport  processes  in  turbulence. 

In  Section  4,  we  report  results  of  Very  Large  Eddy  Simulations  of  turbulent  flow  past 
prototype  propeller  configurations.  In  particular,  we  concentrate  attention  on  the  so-called 
MIT  flapping  foil  experiment  in  which  an  oscillating  NACA16  hydrofoil  is  studied.  This 
example  has  the  advantage  that  it  has  been  widely  studied  previously  using  alternative 
turbulence  models  and  has  an  extensive  and  well  documented  database  associated  with  it. 
The  flapping  foil  experiment  has  been  a  difficult  one  for  previous  modeling  efforts  using 
traditional  two-equation  turbulence  models,  in  the  sense  that  it  has  been  necessary  to 
fix  a  number  of  ad  hoc  modeling  parameters  to  achieve  agreement  with  even  some  of  the 
experimental  data.  In  this  report  we  show  that  the  RNG-based  two-equation  transport 
model  does  an  excellent  job  of  predicting  such  key  turbulent  flow  characteristics  as  skin 
friction  coefficient  Cj  without  adjustable  parameters.  These  results  offer  the  prospect  that 
the  RNG  turbulence  model  will  become  a  standard  for  use  in  Naval  Architecture  community. 
To  this  end,  we  have  worked  to  ensure  that  the  VLES  RNG  model  is  available  in  a  well 
supported,  user-friendly,  generally  available  code  environment.  This  will  help  to  facilitate 
the  broad  application  of  our  work  to  a  variety  of  propeller  design  problems. 
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Chapter  2 


Interaction  of  Surface  Waves  with 
Free  Surface  Turbulent  Flow 


Here  we  describe  numerical  simulations  of  incompressible  unsteady  open-channel  flow.  The 
main  goal  of  this  study  was  to  investigate  the  effects  of  mutual  interaction  of  surface 
waves  with  in-depth  turbulence,  especially  the  influence  of  near-surface  turbulence  upon 
dispersion,  propagation,  and  statistics  of  the  waves. 

2.1  Problem  Background 

Free-surface  turbulent  flows  are  of  genuine  interest  for  ocean  technology.  This  problem 
bridges  surface  wave  phenomena  and  fully  developed  hydrodynamic  turbulence.  For  the 
most  part,  previous  studies  of  surface  wave  phenomena  have  concentrated  on  the  behavior 
of  non-dissipative,  although  highly  dispersive  and  nonlinear,  ensembles  of  waves.  In  this 
case,  potential  flow  is  a  good  approximation  since  molecular  viscosity  is  important  only 
within  a  thin  boundary  layer  j which  is  normally  much  smaller  than  all  the 

characteristic  scales  such  as  the  capillary  length  or  system  dimensions.  It  is  usually  assumed 
that,  inside  this  viscous  layer,  large  vorticity  is  generated.  These  vortex  sheets  move  with 
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the  surface  but  play  a  rather  minor  role  in  the  large-scale  surface  dynamics.  Their  role, 
however,  may  be  important  in  small-scale  ocean-air  interaction,  e.g.  short  gravity  waves, 
capillary  waves  and  wind  generated  ripples.  At  the  same  time,  dissipation  and  vorticity 
generation  phenomena  are  essential  within  the  turbulent  fluid. 

Here,  we  concentrate  upon  the  effects  of  the  mutual  influence  of  random  surface  waves 
and  statistical  properties  of  the  near  free  surface  (FS)  region  in  fully  developed  turbulent 
flow.  To  date,  there  is  no  comprehensive  strong  coupling  theory  of  these  phenomena. 
Numerical  simulations  reported  here  were  done  using  spectral  methods  and  capitalize  upon 
the  success  of  early  numerical  experiments  experiments  for  wall-bounded  flows. 

At  the  present  time,  there  exist  rather  extensive  experimental  studies  of  turbulent  open 
channels  (Ueda  et  al.  1977;  Nezu  k  Rodi  1986;  Kirkgoz  1989),  as  well  as  semi-empirical 
theories  and  results  oi  K  —  £  and  Reynolds  stress  modeling  (Gibson  k  Rodi  1989;  Celik 
k  Rodi  1984;  Swean  et  al.  1991).  The  advantage  of  the  DNS  approach  taken  here  is  its 
independence  of  any  ad  hoc  assumptions  often  necessary  for  theory  to  progress  and  its 
greater  flexibility  in  matching  experimental  conditions.  Our  studies  focused  upon  channel 
flow  bounded  with  a  rigid  wall  from  below  and  having  an  open  upper  surface.  Some  DNS 
calculations  of  open  channel  flow  were  already  performed  in  Handler  et  al.  (1991),  Leighton 
et  al.  (1991)  and  Swean  et  al.  (1991).  However,  in  those  early  works  only  the  case  of  zero 
Froude  number  was  considered.  The  numerical  scheme  used  in  our  simulations  coincides  in 
the  case  of  zero  Froude  number  with  that  of  Swean  et  al.  (1991).  In  our  work  we  mostly 
concentrate  on  the  study  of  the  nonzero  Froude  number  case. 

In  the  following  sections,  we  will  formulate  our  approach  to  the  problem,  then  describe 
the  numerical  scheme,  and  present  the  results  of  flow  simulations. 
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2.2  Strategy  of  Full  Navier-Stokes  Simulation  of  the 
Turbulent  Free- Surface  Problem 

We  solve  the  three-dimensional  Navier-Stokes  equation  for  incompressible  flow: 

dtVi  =  [V  X  —  Vj-TT  -1-  l/AVi 

(2.1) 

P  =  0 

OXi 

where  v  is  the  velocity  field,  p  is  the  pressure,  tt  =  pf  p  +  u^/2,  v  is  the  kinematic  viscosity 
and  p  is  the  density  of  fluid.  The  boundary  conditions  are: 

wall 

(2.2) 

niTjj  surface  '^i'^ij\ext 

where  Tij  =  —p^ij  +  I'P  {diVj  -t-  djVi)  are  the  internal  and  external  components  of  the  stress 
tensor  respective  to  the  fluid,  and  ra,  is  the  normal  to  the  free  surface  (and  summation  over 
repeated  indices  is  implied).  The  system  (2.1),  (2.2)  is  closed  using  the  equation  for  the 
free  surface  height  h: 

dth  4"  Vqi  •  dcth  —  Vnormallfs  (2-3) 

Here  an  index  denoted  by  a  Greek  letter  a  corresponds  to  directions  in  the  free  surface 
(the  X  —  y  plane)  and  da  stands  for  the  differentiation  in  that  direction.  If  we  consider  a 
channel,  the  rigid  wall  is  located  at  z  -  -/f/2,  the  upper  boundary  of  the  physical  domain 
is  given  by  2  =  i//2  -|-  h{x,y,t),  and  the  subscript  fs  designates  values  taken  at  the  free 
surface. 

Our  goal  is  to  use  high-resolution  spectral  methods.  When  this  technology  is  applied  to 
a  flow  bounded  between  two  rigid  walls,  the  flow  is  reproduced  by  a  Fourier  spectral  series 
in  the  x  —  y  direction  and  by  a  Chebyshev  polynomial  representation  in  the  2— direction. 
The  latter  concentrates  collocation  points  towards  the  walls,  providing  high  resolution 


of  boundary  layers.  However,  in  an  open  channel  the  upper  surface  is  curved  and  time- 
dependent,  which  complicates  the  use  of  spectral  methods.  One  possible  approach  is  to  give 
up  the  high  resolution  in  the  direction  provided  by  using  the  Chebyshev  polynomials. 
However,  we  believe  that  it  is  essential  to  ensure  proper  resolution  of  near-free-surface 
boundary  layers. 

The  approach  used  here  is  based  upon  linearization  of  the  boundary  conditions  at  the 
free  surface.  Then  the  boundary  conditions  (2.2)  are  adjusted  to  those  imposed  at  the 
unperturbed  surface.  Such  an  approach  retains  the  ability  to  use  of  the  Chebyshev  spectral 
method  in  the  z-direction  to  achieve  high  resolution.  The  boundary  conditions  then  take 
the  form 


dth  da  {hva)  = 

pU  {da'^z  T  ^z'^cc')  —  j (^*^) 

gh  —  (j^^Ah  +  2vdzVz  =  -  — 

P  P 

where  g  is  the  gravitational  constant,  cr*  =  {cr/p)  is  the  surface  tension,  \ext  and  Pext 
are  the  external  stress  and  pressure,  respectively.  All  terms  in  (2.4)  are  evaluated  at 
h{x,y,t)  =  0;  again  the  index  a  means  that  only  horizontal  components  are  used.  The 
only  nonlinear  term  retained  in  (2.4)  is  the  term  da{hva)  responsible  for  the  convection 
of  the  free  surface  by  the  horizontal  velocity.  It  turns  out  that  in  our  formulation  the 
convective  term  is  the  same  order  of  magnitude  as  though  dah  is  small.  Equations 
(2.1)  and  (2.4)  represent  the  final  formulation  of  the  boundary  problem  that  we  use  here 
,  to  describe  open  surface  turbulent  flow.  In  linear  approximation,  this  system  reproduces 
both  the  dispersion  relation  of  surface  waves  and  the  Helmholtz  instability  at  the  surface. 
Interaction  effects  are  accounted  for  by  the  nonlinearity  of  the  Navier-Stokes  equation.  The 
ability  to  describe  effects  of  nonlinear  self-interactions  of  surface  waves,  a  separate  problem 
which  is  outside  the  scope  of  this  project,  is  sacrificed  when  we  use  the  present  approach. 
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2.2.1  Applicability  limits 

First  of  all,  a  linear  analysis  (Landau  &  Lifshitz  1987)  shows  that,  for  each  spatial  harmonic 
with  a  wave-number  k  and  the  amplitude  h,  linear  approximation  of  the  boundary  condi¬ 
tions  is  only  valid  when  h  <  where  and  il{k)  =  +  It  may 

be  seen  later  that  Fourier  spectrum  of  h{k)  is  rather  steep  and  the  condition  kh{k)  <  lyj 
is  automatically  satisfied  whenever  it  is  valid  at  the  smallest  wave-numbers.  This  leads  to 
the  estimate  for  the  flow  ’integral’  Froude  number  (here  and  later  cr*  =  0  is  assumed  ) 

Fi  =  U/{gHY^^  «  10/i?e^/^  (2.5) 

Here  this  integral  Froude  number  and  the  Reynolds  number  Re  =  U H j i>  are  defined  using 
the  centerline  velocity  U  and  the  channel  height.  Deriving  (2.5)  we  assumed  that  h  ~ 
ul^j2g.  We  also  assumed  that  in  the  range  of  flow  Reynolds  numbers  of  interest  here. 
Re  «  10^  -  10^,  the  characteristic  rms  velocity  at  the  surface  can  be  heuristically  estimated 
&.S  U-pms  ^  IQ-W  (Kim,  Moin  &  Moser  1987). 

At  the  same  time  it  would  not  be  satisfactory  to  consider  only  waves  that  are  shorter 
than  the  viscous  sublayer  width  /„  and  whose  amplitude  are  smaller  than  ly.  The  dissipation 
length  of  waves  ly  is  defined  as  the  scale  where  the  frequency  of  wave  is  equal  to  the  wave 
decay  rate  (see  (2.17)  below).  If  that  was  the  case,  our  analysis  would  only  be  of  general 
methodological  interest  and  not  applicable  to  the  description  of  real  phenomena  in  water 
waves.  Indeed,  in  water  at  normal  conditions,  /„  ~  10“^cm,  whereas  the  capillary  length 
d  =  {2(T»/gY^’^  =  10“^cm.  It  is  clear  that  any  length  scales  smaller  than  d  are  not  of 
interest  for  engineering  and  oceanographic  applications  because  they  are  effectively  damped 
by  capillary  effects. 

The  condition  for  the  existence  of  wave  motion  is  ly  <  fi,  resulting  in  the  inequality 

Fi  >  10/Re^/",  (2.6) 

It  is  clear  that,  as  Re  oo  there  exists  a  wide  range  of  Froude  numbers  where  both 
inequalities  (2.5)  and  (2.6)  hold.  It  is  worth  mentioning  that  in  a  typical  flow  of  engineering 
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interest  {H  ^  Im,  U  ^  O.lm/s,  Re  10^,  Fi  10  the  condition  (2.5)  is  satisfied. 
This  is  also  the  case  for  a  typical  laboratory  setup  in  free-surface-flow  studies  of  Ueda  et 
al.  1977. 


2.3  Spectral  Methods  for  Free  Surface  Turbulence 

A  spectral  method,  with  Fourier  series  in  the  streamwise  (x)  and  spanwise  (y)  directions, 
and  Chebyshev  polynomial  expansion  in  the  normal  direction  (z),  is  employed  spatially. 
A  conventional  time-splitting  scheme  is  used  to  separate  nonlinear  and  viscous  time  steps. 
However,  the  solution  for  incompressibility  and  viscosity  is  done  within  the  same  step 
by  solving  a  fourth-order  equation  for  the  normal  velocity  and  a  second-order  equation 
for  the  normal  component  of  vorticity.  Streamwise  and  spanwise  velocity  components 
are  recovered  from  the  incompressibility  condition.  If  we  define  the  Reynolds  number  as 
Re  =  UHlv,  the  convective  term  as  {);  =  [v  x  w],-  and  the  normal  component  of  vorticity 
as  V’  =  dvxjdy  —  dvy/dx,  (2.1)  takes  the  form: 

|vV  =  /.+  ivV.  (2.7) 

=  +  (2.8) 

where 

d 

fy  =  Av^  -  —div{v) 

and 

=  rot^{v). 

Here  the  incompressibility  condition  has  already  been  taken  into  account. 

The  boundary  conditions  at  the  wall  are: 

0 

0  Itua//  0  (2"^) 

The  boundary  conditions  at  the  free  surface  are  more  complicated.  The  conditions  on 
the  tangential  components  of  stresses  (see  (2.4))  can  be  expressed  straightforwardly  in  the 
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form: 


(2.10) 


dx 

q2  (3  3  \ 

I/s  l/«~  '^xz  \ext  '\'~^Tyz  \extj  ~  ^div  (^’H) 

where  the  subscript  a  again  indicates  horizontal  components.  The  third  boundary  condition 
on  Vz  is  obtained  through  the  relation  on  pressure  at  the  surface.  In  our  numerical  scheme, 
the  pressure  computation  is  not  required  for  time  advancement.  Should  data  on  pressure 
be  needed  for  turbulent  statistics,  it  can  be  calculated  by  two  different,  but  essentially 
equivalent  ways.  One  can  either  use  the  equation  for  the  normal  component  of  velocity 
with  the  wall  pressure  values  determined  from  an  equation  for  daVa,  or  one  can  use  the 
equation  for  daVa  with  the  pressure  corresponding  to  zero  wave-numbers  determined  from 
the  u^-equation.  Therefore,  it  is  convenient  to  impose  the  pressure  boundary  condition  in 
the  following  form: 


3 

—tp  |/,=  Re 


‘yz 


_l_  I 

ext  o  Trz  \ext  I  —  ^rot 

dy  J 


3  1 

■^dzVz  =  A„7r  -  daVa  +  ^  A  (d^Vz) 


2 

TT  =  gh-  a*Ah  -f-  -^OzVz  +  -77-  Pext 
lie  L 


(2.12) 


Here  the  equation  for  dvzidz  is  obtained  from  the  equation  for  daVa- 

Time  stepping  is  carried  with  a  semi-implicit  scheme  involving  a  Crank-Nicolson  scheme 
for  the  viscous  terms  and  Adams-Bashforth  scheme  for  the  nonlinear  terms.  For  the  time 
advancement  of  surface  height  we  use  a  Crank-Nicolson  scheme  for  the  normal  velocity 
component  and  an  Adams-Bashforth  scheme  for  the  convective  term  (Canute  et  al.  1987). 
Equation  (2.8)  then  reduces  to 

('  -  =  f  +  ('  +  2!^^ 

with  the  boundary  conditions  (2.9),  (2.10).  Equation  (2.13)  is  solved  by  the  Chebyshev-tau 
method  for  each  horizontal  Fourier  component  of  the  vorticity  field. 
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The  fourth-order  equation  (2.7)  is  solved  by  splitting  it  into  two  second-order  equations 


1  - 


6t 

2Re 


U 


n+1  _ 


a 

2 


(3/,”-/r')  +  (i  + 


(2.14) 


=  4, 


n+1 


with  the  boundary  conditions  (2.9),  (2.11)  and  the  pressure  boundary  condition  (2.12) 
which  has  the  following  form: 

Here  we  use  the  following  notations: 


f^n+l/2  ^  1/2(A”+1  +r), 

g*=g  +  <T*{kl  -I-  k])  and 


g  —  6t  ^  doiVct  2  )  • 


The  time  advancement  scheme  for  surface  height  is 


8t 

2 


h  = 


d 

dXa 


(UQ./i(XQ.  )  ) 


(2.16) 


The  system  of  coupled  equations  (2.14),  (2.16)  and  (2.9),  (2.11),  (2.15)  is  solved  by  the 
Chebyshev-tau  method,  in  which  the  four  boundary  conditions  are  satisfied  by  means  of  a 
Green’s  function  technique.  This  system  requires  three  Green  functions  Gi-  The  first  two 
of  them  solve  the  equations 

(i-+vdft,  =  0 


2Re 


V^Gi,2  —  Gi.2 
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with  the  boundary  conditions 


and 


AGi  l/s=  Gi  I/s—  Gi  Iwaii—  0; 
AGi  |tya//=  1; 

AG2  —  G2 

G2  |/s=  1; 


52^2 

dz^ 


|/s=  -k‘^G2  I/s  ■ 


In  order  to  take  into  account  the  ’tau  correction’,  i.e.  numerical  errors  due  to  the  finite 
number  of  retained  Chebyshev  polynomials,  it  is  necessary  to  introduce  a  third  Green 
function  G3  which  satisfies  the  equation 


AG3  =  0 


with  the  boundary  condition 


G3  \wall  —  0) 

G3  I/s  =  1- 

The  third  Green  function  G3  is  used  to  satisfy  equations  at  the  boundary  points  and  the 
boundary  conditions  at  the  same  time. 

To  avoid  aliasing  errors  involved  in  computing  the  nonlinear  terms  pseudo-spectrally, 
the  usual  2/3  dealiasing  procedure  is  used. 

The  computations  reported  here  were  performed  on 

The  accuracy  of  the  numerical  code  was  examined  in  two  different  ways.  First  we 
implemented  the  numerical  scheme  for  closed  channel  flow  and  analyzed  the  time  evolution 
off  small-amplitude  Orr-Sommerfeld  modes.  In  all  the  tests  performed  both  linear  decay 
and  growth  rates  were  predicted  with  errors  of  less  then  10~^%.  In  the  second  set  of  tests, 
we  analyzed  the  decay  in  time  of  an  infinitesimally  small  wave  with  given  wave  number  k 
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at  the  free  surface.  We  measured  the  frequency  0,(k)  and  the  decay  rate  r(A;)  of  the  wave 
and  compared  them  with  the  theoretical  prediction  for  small  amplitude  deep  water  waves 
(Landau  &;  Lifshitz  1987) 

n(k)  =  {gk  +  Y{k)  =  2uk\  (2.17) 

Equation  (2.17)  is  valid  when  F  <<  fl,  which  is  always  the  case  for  our  range  of  parameters. 
It  was  shown  that  the  frequency  Q.{k)  is  reproduced  by  the  numerical  scheme  very  accurately 
(usually  in  the  range  of  10“^  -  10"^%),  provided  the  time  step  St  is  small  enough  Q.{k)St  < 
0.2.  Our  numerical  time  advancement  scheme  has  second  order  accuracy  in  time.  Since 
r  <<  it  is  necessary  to  decrease  St  even  further  if  we  want  to  obtain  the  correct  value 
for  the  decay  rate  of  waves.  It  may  be  shown  that  significant  numerical  dissipation  may  be 
avoided  if  StQ.{k)  <  {T{k)l^{k)Yl'^  {Ivkf/'^.  Under  this  restriction,  which  is  generally 

not  so  burdensome,  we  obtained  a  decay  rate  that  differs  from  the  theoretical  one  by  only 
10"^%.  In  all  our  runs,  we  chose  the  time  step  small  enough  to  correctly  reproduce  the 
dynamics  of  the  surface  waves. 

2.4  Computational  Data 

We  consider  a  setup  where  turbulence  is  generated  through  an  externally  imposed  pressure 
gradient  in  the  x  direction  which  supplies  the  centerline  mean  velocity  U .  The  initial  ve¬ 
locity  fields  were  obtained  from  the  standard  Orr-Sommerfeld  instability  modes  for  channel 
flow.  The  flux  in  the  x  direction,  i.e.  the  bulk  mean  velocity 

1  /•!  ,2 

^  2  i-i 

was  fixed  by  adjusting  the  mean  pressure  gradient.  The  viscosity  is  chosen  to  be  —  1/3000 
and  the  maximum  mean  velocity  Umax  ~  0.765.  The  Reynolds  number  defined  through  the 
maximum  mean  velocity  is  Re  ^  4560.  The  wall-shear  velocity  is 

u.  =  «  0.041 
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so  that  i?e*  =  {u^H)/^  ~  250. 

Under  these  conditions,  the  Kolmogorov  dissipation  wave-number  estimated  via  the 
energy  dissipation  rate  at  the  centerline,  ~  5  x  is  kd  =  Ri  35  and  is 

smaller  than  the  maximum  wave-number  in  the  x  —  y  plane.  This  means  that  turbulence 
in  the  bulk  of  the  channel  and  near  the  rigid  wall  is  well  resolved.  We  believe  that  the 
near-free- surface  region  is  also  resolved  properly.  Indeed,  the  Froude  number  Fi  in  our  runs 
varied  from  0.1  to  0.6,  and  the  corresponding  boundary  layer  was  thick  enough  to  contain 
several  Chebyshev  collocation  points,  (although  /„  is  less  than  the  Kolmogorov  dissipation 
scale).  Both  criteria  (2.5)  and  (2.6)  hold  in  this  case  and  the  linearization  of  the  boundary 
conditions  is  well  based. 

2.4.1  Turbulence  statistics 

The  mean  velocity  profile  is  shown  in  Figure  1.  In  Figure  1,  the  velocity  is  normalized  by  the 
wall-shear  velocity  n.  and  the  distance  from  the  rigid  wall  z+  is  measured  in  wall  units  = 
As  expected  the  mean  velocity  profile  obeys  the  law  of  the  wall  =  2.blnz^  +  5.5 
starting  from  approximately  z+  =  30.  Typically,  the  mean  velocity  is  slightly  larger  than 
one  would  expect  from  the  law  of  the  wall  just  near  the  free  surface,  in  accordance  with  the 
data  Nezu  &  Rodi  (1986).  Such  deviations  are  likely  to  result  from  the  effects  of  the  free 
surface.  The  data  plotted  in  Figure  1  corresponds  to  a  typical  run  with  Froude  number 
Fi  «  0.55. 

Turbulence  intensities  and  Reynolds  shear  stress  normalized  by  the  wall-shear  veloc¬ 
ity  are  shown  in  Figures  2  and  3,  respectively.  The  data  on  rms  vorticity  is  shown  in 
Figure  4.  Near  the  rigid  wall  and  close  to  the  centerline,  the  profiles  are  in  good  agree¬ 
ment  with  the  DNS  data  on  rigid-wall  channel  flow  (Kim,  Moin  &  Moser  1987)  and  in 
reasonable  agreement  with  experimental  data  (Kreplin  &  Eckelmann  1979).  For  example, 
the  intensity  of  streamwise  turbulent  fluctuations  peaks  at  z+  13.  We  have  compared 
our  measurements  for  the  skin  friction  coefficient  Cj  =  T^aiil\U^  with  the  experimental 
correlations  proposed  by  Dean  (1978).  The  friction  coefficient  in  our  case  is  found  to  be 
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Cj  ~  7.6  X  10  This  agrees  well  with  Dean’s  formula  Cf  ~  0.073jRe^°'^®  ~  7.7  x  10  ^ 
where  Re^  =  2HUmli'  ~  8000.  The  ratio  UmaxIUm  is  1.148,  also  in  a  good  agreement 
with  Dean’s  data  UmaxIUm  =  1.28i?e;;°  “^^®  =  1.153.  The  comparison  with  available  data 
on  turbulence  intensities  in  open  channel  flow  (Nezu  &  Rodi  1986)  also  shows  good  agree¬ 
ment.  The  streamwise  rms  velocity  at  the  free  surface  in  wall  units  is  roughly  1  in  our 
study  and  0.8  —  1.2  in  the  experiments  and  the  spanwise  rms  velocity  is  roughly  0.7  and 
0.65,  respectively. 

An  important  outcome  of  our  study  is  that  within  the  feasible  range  of  Froude  numbers 
(2.5),  typical  mean  flow  profiles  do  not  change  significantly  with  the  flow  Froude  number. 
Turbulence  intensities  and  rms  vorticity  data  are  also  weakly  sensitive  to  the  Froude  num¬ 
ber.  We  may  conclude  that  for  small  Froude  numbers  it  is  reasonable  to  consider  even 
the  limiting  case  Fi  -  0.  Such  flow  in  which  the  gravitation  constant  g  is  infinite  can 
be  visualized  as  corresponding  to  a  half-channel  flow  with  two  rigid  walls,  where  the  in- 
viscid,  no-stress  boundary  conditions  are  imposed  at  the  wall  corresponding  to  the  ’free 
surface’  (Hunt  &:  Graham  1978).  The  presence  of  such  a  no-stress  wall  is  the  dominant 
effect  at  small  Froude  numbers.  In  this  sense,  the  free  surface  responds  nearly  linearly  to 
the  bulk  turbulence  and  does  not  lead  to  significant  changes  in  the  fluid  flow  near  the  free 
surface  compared  with  the  case  Fi  =  0.  The  dependencies  on  the  flow  Froude  number 
are  only  significant  in  the  near-free-surface  region  of  the  size  =  {2v where  fi/i  is 
some  characteristic  frequency  of  turbulence  at  the  free  surface.  These  effects  are  discussed 
further  in  the  next  Section. 

The  most  important  influence  of  the  free  surface  is  that  the  fluid  motion  near  the  free 
surface  is  quasi  two-dimensional.  The  fluctuations  of  the  velocity  normal  to  the  surface 
are  substantially  smaller  than  the  ones  in  horizontal  directions.  The  only  component  of 
vorticity  that  is  large  near  the  free  surface  is  that  normal  to  the  free  surface.  We  may  expect 
that  the  most  important  excitations  of  the  free  surface  come  from  vortex  tubes  attached 
to  the  free  surface  at  one  end,  produced  by  small  scale  eddies  of  nearly  two-dimensional 
turbulence.  In  this  case,  an  inverse  energy  cascade  may  take  place  that  may  facilitate  the 
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formation  of  large  vortex  tubes  near  the  free  surface. 


Turbulent  eddy  viscosity  in  the  near  free  surface  region 

Using  data  on  the  Reynolds  stress  and  mean  velocities,  we  obtained  distributions  of  eddy 
viscosity  in  open  surface  flow.  A  typical  profile  of  eddy  viscosity  is  plotted  in  Figure  5  at 
Fi  =  0.55. 

First,  we  compared  our  data  with  the  empirical  formula 


^eddy 


z  ,  z  . 

0A-(l  -  -). 


(2.18) 


A  comparison  of  (2.18)  with  experimental  data  for  open  channel  flow  was  made  by  Ueda 
et  al.  (1977).  Eq.  (2.18)  implies  that  the  law  of  the  wall  is  valid  throughout  the  channel, 
which  is  certainly  not  true  at  least  near  the  rigid  wall.  From  Fig.5,  we  see  that  the  measured 
distribution  is  rather  well  described  by  (2.18)  near  the  free  surface. 

We  have  analyzed  the  dependence  of  eddy  viscosity  upon  the  turbulent  kinetic  energy 
K  and  the  energy  dissipation  rate  £.  Our  recent  developments  in  renormalization  group 
turbulence  modeling  suggest  that  the  mean  rate-of-strain  dependence  of  eddy  viscosity  is 
also  important  (Yakhot  et  al.  1992).  We  found  the  following  expression  to  be  in  good 
agreement  with  the  present  data: 


2 

5  1  +  {tj/ay 


(2.19) 


Using  turbulent  characteristics  K  and  5  measured  from  the  simulations,  we  find  that  the 
best  choice  of  numerical  constants  is  =  0.09,  a  =  3.0.  The  dimensionless  parameter  rj 
is  the  ratio  of  turbulent  to  mean  strain  time  scale.  Near  the  rigid  wall  7/  is  large  (~  25); 


in  the  log-layer  of  the  channel  rj  is  approximately  constant  (between  3  and  4);  near  the 
free  surface  rj  approaches  zero.  As  the  results  plotted  in  Fig.  5  suggest,  the  interpolation 
formula  (2.19)  holds  rather  well  near  the  rigid  wall.  Close  to  the  free  surface,  the  eddy 
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viscosity  given  by  (2.19)  differs  significantly  from  the  measured  data,  due  to  the  decrease 
of  energy  dissipation  rate  caused  by  the  zero  stress  boundary  conditions. 


2.4.2  Details  of  near- free- surface  behavior 


Our  results  show  that  the  viscous  term  in  the  boundary  condition  (2.12)  is  numerically 
small  and  prms  ~  ghrms-,  where  Kras  -<  is  the  rms  surface  height.  On  the  other 

hand,  the  rms  pressure  at  the  free  surface  prms  ?«<  u?  >  /2  is  approximately  independent 
of  g  (see  Fig.  6).  Therefore,  the  following  approximate  relation  is  well  satisfied: 

A™.  «  (2.20) 
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where  K js  is  the  mean  kinetic  energy  at  the  surface. 

Another  observation  is  that  the  root  mean-square  fluctuations  of  the  streamwise  and 
spanwise  velocity  components,  as  well  as  the  normal  component  of  vorticity  are  practi¬ 
cally  independent  of  g,  whereas  |rm5~  ^/g-  The  rms  values  of  streamwise  and  spanwise 
components  of  vorticity  are  also  inversely  proportional  to  g.  The  scaling  of  with  g  and 
(2.20)  suggest  that  there  exists  a  characteristic  frequency  of  surface  height  fluctuations 
\rms^  '^z  |rms  /^rms  which  is  approximately  constant  in  the  range  of  Froude  numbers 
considered.  This  frequency  relates  to  the  characteristic  turnover  times  of  the  largest  turbu¬ 
lent  eddies.  Within  the  viscous  sublayer,  the  turbulent  kinetic  energy  K  was  observed  to 
be  practically  independent  of  the  distance  from  the  free  surface.  Therefore,  the  Neumann 
boundary  condition  dKjdz  ==  0  seems  to  be  relevant  for  modeling  free-surface  phenomena. 

Most  of  the  averaged  quantities  or  their  derivatives  vary  smoothly  in  the  middle  of  the 
channel  but  vary  nearly  exponentially  near  the  free  surface.  To  analyze  their  behavior,  it 
is  convenient  to  define  an  operator  T  of  the  following  form: 


/(^/»)  ~  /(°°) 

/(O)  -  /(oo)  ■ 


(2.21) 


The  operator  T  acts  on  functions  f{zfs),  where  zjs  designates  the  distance  from  free  sur¬ 
face.  According  to  the  definition,  =  1  for  any  /.  For  functions  which  fall  off 
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exponentially  away  from  the  free  surface,  ~  ^xp(—ZfsfX).  In  Figure  7,  we  show 

the  near-free-surface  behavior  of  different  flow  characteristics  analyzed  by  means  of  the 
limiting  operator  J-.  Plotted  in  Figure  7  are  the  mean  rate  of  strain  <  dvxjdz  >,  the 
energy  dissipation  rate  the  horizontal  components  of  rms  vorticity  as  well  as  the 

z-  derivatives  of  the  Reynolds  stress  <  d{vxV:,)ldz  >.  Most  of  these  characteristics  change 
exponentially  near  the  free  surface  as  ~  Ci  —  C2Cxp[—ZfsfX)  with  approximately  the  same 
exponent  A  0.0877'.  If  we  identify  A  with  a  depth  of  some  skin-layer,  a  characteristic 
excitation  frequency  at  the  free  surface  can  be  estimated  as  fio  ~  2ufX^.  In  our  setup, 
such  estimation  gives  fio  ~  0.42.  It  is  shown  later  that  this  frequency  nearly  matches  the 
directly  measured  characteristic  frequency  of  surface  wave  motion. 

As  may  also  be  seen  from  Figure  7,  the  ratio  of  the  turbulent  kinetic  energy  K{zfs)l K{0) 
and  the  ratio  of  rms  normal  vorticity  fluctuations  ca™®(2/s)/a;™®(0)  are  approximately 
constant. 

2.4.3  Surface  waves  on  top  of  turbulent  flow 

Dispersion  relations. 

Let  us  define  two  correlation  functions  of  some  variable  q{x,t) 

Uq{k)  =  J  <  q{x,t)q{0,t)  >  exp{ikx)dx  (2.22) 

Sq{k,  n)  =  j  <  q{x,  t)q{0, 0)  >  exp{ikx  -t-  m)dxdt  (2.23) 

so  that  Uq  =  f{d£lf2‘7r)Sq.  Here  k  stands  for  the  wave-vector  in  the  horizontal  plane  and 
all  quantities  are  evaluated  at  the  free  surface.  The  correlation  functions  Uq  are  directly 
measured  in  our  numerical  experiments,  whereas  the  evaluation  of  the  functions  Sq  requires 
rather  expensive  time  averaging.  In  our  case,  a  linearized  version  of  (2.3)  yields: 

=  (2.24) 

and  for  the  equal-time  spectral  densities; 

W  =  (2.25) 
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In  the  low-k  region  (A:  <<  !//„,),  the  singularities  of  correlation  function  Sh  are  dominated 
by  the  oscillatory  part  of  dispersion  relation.  Using  Equation  (2.17),  we  obtain 

UUk)  =  C(k)U..m,  =  (2.26) 

At  this  point,  it  is  important  to  mention  that  the  spectra  Uh  and  Uz  themselves  may  be 
far  from  those  corresponding  to  the  linear  surface  waves  problem.  Also  the  relation  (2.26) 
does  not  require  potentiality  of  the  surface  flow.  The  only  assumption  that  was  used  in 
the  derivation  of  (2.26)  is  that  the  Green  function  for  the  equation  for  normal  velocity  has 
a  pole  at  the  characteristic  frequency  given  by  (2.17),  whereas  the  dissipative  term  in  the 
dispersion  relation  is  negligible.  Such  an  assumption  may  become  invalid  due  to  the  effects 
of  turbulence.  Also,  the  effects  of  the  mean  flow  in  our  case  can,  in  principle,  complicate 
the  dispersion  relation  structure.  The  ratio 

may  be  called  a  characteristic  frequency  of  surface  height  excitations  at  wave  number  k. 
Measurements  of  the  spectra  (2.27)  can  serve  as  a  diagnostic  tool  to  analyze  the  properties 
of  near-free-surface  turbulence. 

In  Figure  8,  we  show  the  two-dimensional  spectrum  of  the  logarithm  of  characteristic 
wave  frequencies.  The  shape  of  the  spectrum  is  practically  isotropic  in  the  kxky—  plane 
so  that  the  dispersion  relation  can  be  well  represented  by  the  ’one-dimensional’  spectrum 
shown  in  Figure  9.  The  one-dimensional  spectrum  is  obtained  from  the  isotropic  two- 
dimensional  spectrum  by  summing  all  the  Uq{k)  with  A:  on  a  circle  of  a  given  radius  [A:]. 
When  k/kd  >  1,  the  measured  dependence  Clh{k)  is  well  described  by  (2.17).  At  smaller 
wave  numbers,  the  function  ilh{k)  deviates  significantly  from  (2.17)  and  reaches  an  ap¬ 
proximately  constant  value  flo-  This  behavior  of  the  characteristic  frequency  reflects  the 
influence  of  ’bottom-generated’  turbulence  upon  the  surface  wave  behavior.  The  wave  spec¬ 
tra  in  the  low  k  region  k/kd  ^  1  forced  by  the  turbulence.  The  frequency  flo  depends 
on  the  flow  Reynolds  number  only  and  characterizes  the  frequency  of  the  attached  vortex 
tubes.  In  the  turbulence  dissipation  range,  nonlinear  effects  are  small  and  surface  waves 
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can  be  considered  as  freely  evolving.  It  should  be  noted  that  the  numerical  value  of  the 
limiting  frequency  Oq  =  ^{k  — 0)  «  0.5  is  close  to  the  estimate  obtained  in  Section  4.2. 

Spectral  exponents. 

The  spectra  of  surface  height  h{k),  and  velocity  components  Vx^y{k),Vz{k)  are  also  isotropic 
in  wave-number  space.  The  one  dimensional  spectra  are  plotted  in  Figure  10.  Spectral 
behavior  of  the  surface  height  correlation  function  differs  significantly  from  the  Phillips’ 
surface  wave  spectral  laws  1/k^  (Phillips  1977)  or  (Phillipsl  1985).  In  fact,  the  spec¬ 

tral  index  is  between  one  and  two.  This  fact  is  not  very  surprising,  because  we  are  mostly 
concerned  with  spectra  of  small  surface  waves  forced  by  the  bulk  turbulence.  Linearization 
of  the  free  surface  boundary  definitely  contribute  to  this  discrepancy.  Surface  waves  in  our 
problem  play  the  role  of  a  diagnostic  tool  that  enable  us  to  measure  spectral  characteristic 
of  bulk  turbulence.  On  the  other  hand,  Phillips’  spectrum  may  be  observed  in  the  so-called 
wave  saturation  range  for  waves  in  the  ocean,  where  nonlinear  wave  interaction  is  essential. 
Our  data  are  in  reasonable  agreement  with  the  low  Froude  number  data  of  Brumley  & 
Jirka  (1987)  where  the  spectra  of  free-surface  velocity  fluctuations  were  measured  for  the 
case  of  grid  generated  turbulence  with  a  grid  at  the  bottom  of  the  tank. 

For  all  Froude  numbers  that  we  considered,  the  characteristic  surface  frequencies  Ctk 
were  lower  than  the  free-surface  frequencies  (2.17).  When  the  gravitational  constant  g 
decreases,  these  frequencies  match  and  resonant  excitations  of  the  surface  waves  occur. 
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Figure  2.1:  Mean-velocity  profile,  2:+  =  0  corresponds  to  the  rigid  wall;  the  dashed  line  is 
the  law  of  the  wall  u+  =  5.5  -t-  2.5lnz^]  Umaxlu*  =  18.6,  UmaxIUm  =  1-148 
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Figure  2.2:  Root-mean-square  velocity  fluctuations  normalized  by  the  wall  shear  velocity 
X,  u’’’”®;  □,  A,  (a)  in  global  coordinates;  (b)  in  wall  coordinates. 
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Figure  2.3:  Reynolds  stress  -  <  >:  (a)  in  global  coordinates;  (b)  in  wall  coordinates. 

Correlation  coefficient  —  <  (c)  in  global  coordinates. 


Figure  2.7:  Near-free-surface  behavior  of  flow  characteristics.  The  limiting  operator  J~  is 
applied  to:  x,  A,  <  d{v:,v^)ldz  >;  □,  <  dv^jdz  >;  o,  The  ratios  of  K{zfs)IK{^) 

and  i^z f s) ! are  shown  by  *  and  •,  respectively.  Solid  line  is  with  A  ~ 

O.OSH.  Zfs  is  the  distance  from  the  free  surface. 


Figure  2.10:  One  dimensional  spectra  of  surface  height  correlator  Uh{k)lh1^^  (solid  line); 
the  normal  velocity  at  the  surface  correlator  Uv^k)!  <  >l/s  (dashed  line);  the  tangential 
velocity  at  the  surface  correlator  ^^{k)/  <  v\  y  >|/s  (dot-dashed  line)  as  the  function 

of  kjk^.  Dot  line  is  Ijk^l^. 


Chapter  3 


Computational  Methods 


3.1  Variational  formulation  of  spectral  elements  method 
for  flows  with  moving  boundaries 

Here  we  outline  the  technical  details  involved  in  implementation  of  the  free  moving  bound¬ 
ary  problem.  We  consider  simulation  of  incompressible  flows  involving  free  surfaces  which 
can  be  moving  in  time.  The  equations  of  motion  are  the  incompressible  Navier-Stokes 
equations: 


—  -f  (V  •  v)v  =  V  ■  cr  =  — Vp  +  nV^v  (3-1) 

V  •  V  =  0  (3.2) 

In  order  to  develop  spatial  discretization,  let  us  derive  the  variational  form  of  the  Navier- 
Stokes  equations  when  the  whole  mesh  (or  part  of  the  mesh)  is  moving  with  a  velocity 
w.  Consider  the  motion  of  fluid  in  a  time-dependent  control  volume  Q{t)  moving  with  the 
velocity  w.  The  mesh  velocity  w  is  chosen  so  that  it  minimizes  the  deformation  of  the 
domain  fl.  According  to  Reynolds  transport  theorem: 

4-  /  (t/>  *  v)  dcu  =  /  ^  {ip  *  v)  doj  +  [  (n-w)  {ip,*  v)  ds  (3.3) 

dt  jQ(t)  Jn(t)  at  JdQ(t) 
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Here,  ip  is  a  vector  test  function  and  *  denotes  collocation  products  where  no  summation 
is  assumed.  This  way,  one  obtains  from  the  first  right  hand  side  term  that 

9  f  I  dv  dip  .  . 

The  first  term  on  the  right  side  of  equation  (3.4)  can  be  substituted  from  equation  (3.2), 
whereas  the  second  term  on  the  right  side  dip/dt  is  only  due  to  the  motion  of  the  domain 
with  velocity  w  and  can  therefore  be  replaced  by 

^  =  ^  +  w-VV>  =  0  (3.5) 

dt  dt 

This  way  we  obtain  from  (3.3),  (3.5),  and  the  divergence  theorem  that 

—  /  (ip*v)  dw  —  f  Ip  *  -^du  —  f  V  *  (w  •  '\/ip)duj  +  f  V  •  {w{ip  *  v))  du  (3.6) 

dt  jQit)  Mt)  dt  JQ{t)  Jm 

and  by  expanding  the  second  term  on  the  right  and  substituting  from  equation 

—  /  (ip*v)(ko  =  f  Ip  •  (a  -  vv)  du)  -  f  v  *  {-w -Vip)  duj  (3.7) 
dt  J^(ty  JQit) 

+  /  Ip  *  CV  •  wv)  du>  /  V  *  (w  •  S/ip)  du) 

jQ(t)  JQ{t) 

and  by  integration  by  parts  of  the  first  right  hand  side  term, 

—  /  (ib*v)du}=  f  Ip  *  V  ■  (wv  —  w)  dui  —  [  (j-Vipdio  +  f  {n-a)ip^ds  (3.8) 

dt  jQity  JQ(t)  JdQ(t) 

which  is  the  variational  form  of  the  flow  equations  for  a  domain  Q,(t).  The  spatial  dis¬ 
cretization  then  proceeds  by  defining  the  test  functions  ip  to  be  the  Legendre  Lagrangian 
interpolants  used  for  Legendre  spectral  elements,  and  by  obtaining  the  discretized  set  of 
equations. 
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3.1.1  RNG  large-eddy  simulations 

Using  this  approach,  we  have  simulated  a  number  of  massively  separated  flows  in  curved 
geometries.  An  efficient  spectral  element- Fourier  code  was  developed  and  the  RNG  subgrid 
model  was  incorporated  into  it  in  order  to  provide  a  platform  for  large  eddy  simulations. 
Flow  past  a  spherical  surface:  Validation  tests  were  performed  in  the  low  Reynolds 
number  regime  where  the  flow  exhibits  the  large  scale  unsteadiness  persisting  in  the  high 
Reynolds  number  regime,  and  the  results  compare  very  favorably  with  experimental  results 
{Re  <  1,000).  In  the  high  Reynolds  number  regime,  large  eddy  simulation  was  performed 
at  Re  =  20, 000.  This  flow  is  characterized  by  a  large  scale  vortex  street  and  a  complex 
small  scale  flow  structure  in  the  near  wake.  The  time  behavior  of  the  drug  coefficient 
was  highly  unsteady  with  the  characteristic  frequency  of  approximately  twenty  times  the 
Strouhal  frequency  in  agreement  with  the  available  experimental  data. 

3.2  Implementation  of  Transport  models 

We  have  implemented  the  renormalization  group  theory  K  —  e  transport  model  on  a  finite 
volume  platform  which  handles  arbitrarily  complex  geometries  using  curvilinear  coordinate 
systems.  The  main  advantage  of  using  a  finite  volume  discretization  is  that,  being  a  low 
order  method,  it  is  very  efficient,  in  terms  of  both  speed  and  memory.  In  addition,  since 
the  main  purpose  of  turbulence  transport  modeling  is  to  solve  for  an  average  field  where 
all  small  scale  fluctuations  in  space  and  time  are  smoothed  out,  the  optimal  dispersive 
and  diffusive  properties  of  spectral  elements  (or  other  spectral-type  global  methods)  are 
not  really  needed.  On  the  other  hand,  because  of  the  global  character  of  the  spectral-type 
methods  which  couple  many  points  together,  it  is  extremely  difficult  to  implement  special 
treatment  close  to  walls  when  using  low  Reynolds  number  models  which  are  essential  for 
the  non-equilibrium  regions.  Finite  volume  methods  allow  for  local  use  of  wall  layers  and 
direct  control  of  grid  size,  necessary  for  the  low  Re  number  wall  models. 

Even  in  the  case  of  large  eddy  simulation  (LES),  where  no  separate  equations  are  solved 
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for  the  turbulent  kinetic  energy  and  energy  dissipation  rate,  there  is  certain  ambiguity  in 
the  definition  of  the  resolvable  flow  length  scale  (the  cut-off  of  filter  function)  in  spectral- 
type  methods.  This  is  a  direct  consequence  of  the  global  character  of  spectral  methods  and 
the  absence  of  a  well  defined  local  grid  size,  in  contrast  with  the  finite  volume  methods, 
where  the  grid  size  is  obvious. 

3.2.1  Formulation  of  the  model 


The  incompressible  Reynolds-Averaged  Navier-Stokes  (RANS),  closed  by  the  RNG  JC  —  £ 
model  is  given  by 


dUi 

dt 


+  UjV.U, 


ViUi  =  0 


^  -h  t/jVj/C  =  UeddyS'^  -  £  + 


f  -b  UjVjS  =  C,,  ^UeddyS^  -C,,^-R  +  VavT^e 


where  vj  =  I'eddy  +  ^moi  and  the  rate  of  strain  term  R  is  given  by 


R  — 


V 


dui  dui 
dxi  dxj 


which  is  expressed  in  RN  G  }C  —  £  model  as 

1  +  ^7/3  fc 

where  t]  =  S]C/£  and  =  2SijSij  is  the  magnitude  of  the  rate-of-strain. 

These  equations  are  incorporated  in  a  code  which  solves  the  governing  partial  differential 
equations  for  the  conservation  of  mass  and  momentum,  as  well  as  the  transport  of  turbulent 
kinetic  energy  and  the  rate  of  dissipation.  A  general  form  of  the  above  PDEs  can  be  written 
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in  Cartesian  tensor  notation  as 


_a 

dt 


(«)  +  £:(«.?) 


dxi  dxi 


+  Sq 


where  q  is  the  conserved  quantity.  The  equations  are  discretized  by  integration  over  the 
grid  cells  which  results  in  a  set  of  algebraic  equations  of  the  following  general  form: 


QP  ~  ^p)  — 

i  i 

where  the  summation  is  over  the  neighboring  cells.  The  A—  coefficients  contain  contribu¬ 
tions  from  the  convective  and  diffusive  fluxes  whereas  Sc  and  Sp  are  the  components  of 
the  linearized  source  term,  5,  =  Sc+Spqp.  A  second  order  power  law  difference  scheme  is 
used  for  the  interpolation  between  grid  points  and  calculation  of  the  derivatives.  The  set 
of  corresponding  algebraic  equations  is  solved  by  Patankar  semi-implicit  iterative  schemes 
using  relaxation.  The  sequence  of  iteration  steps  is  as  follows: 

1.  The  momentum  equations  are  solved  using  current  values  for  pressure,  in  order  to 
update  the  velocity  field. 

2.  The  Poisson  equation  is  derived  from  the  continuity  equation  and  the  linearized 
momentum  equations.  This  “pressure  correction”  equation  is  then  solved  to  obtain  the 
necessary  corrections  to  the  pressure  and  velocity  fields  such  that  the  continuity  is  achieved. 

3.  The  K  and  e  equations  are  solved  using  the  updated  velocity  and  the  current  value 
of  eddy  viscosity.  Then  the  eddy  viscosity  is  updated. 

4.  Check  for  convergence.  Convergence  is  reached  when  the  sum  of  normalized  residuals 
falls  below  certain  value.  For  a  general  variable  q,  the  normalized  residual  is  defined  as: 


^  _  Y^nodesP  lAcgg  +  Awqw  +  AjvgAf  +  Asqs  +  Sc  —  Apqp\ 

'YnodesP  lAp^pj 

where  Ap,  Aw,  An^As  are  the  coefficients  which  combine  convection  and  diffusion  through 
the  control  volume  surrounding  given  point  P.  Sc  and  Sp  are  the  components  of  the 
linearized  source  term.  The  summation  here  extends  over  all  of  the  computational  points. 
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Chapter  4 


RNG  Transport  Modeling  of 
Hydrofoil  Flows 


In  our  simulations  of  complex  hydrofoil  flow  using  the  RNG  transport  model  we  have 
observed  excellent  agreement  with  the  experimentally  measured  data,  including  the  skin 
friction  and  surface  pressure  coefficients,  the  location  of  separation  region  near  the  trailing 
edge  of  the  testfoil,  the  wake  velocity  profiles  and  the  boundary  layer  profiles  on  the  foils.  At 
the  same  time,  the  RNG  transport  model  has  demonstrated  significant  superiority  over  the 
standard  standard  K  —  t  model  which  does  equally  good  job  on  the  mean  flow  prediction  but 
overestimates  the  friction  coefficient  by  a  factor  of  2.  Here  we  will  describe  the  experimental 
setup  which  was  fully  matched  in  our  simulations,  outline  the  numerical  techniques,  and 
analyze  the  computational  data. 


4.1  Outline  of  the  Benchmark  Experimental  Setup 


Experiments  have  been  performed  in  the  MIT  Variable  Pressure  Water  Tunnel  to  study  the 
flow  over  an  isolated  two  dimensional  hydrofoil  subject  to  steady  loading  and  to  vertical 
gust  loading  at  high  reduced  frequency.  A  modified  NACA16  hydrofoil  with  chord  length  of 
18  inches  and  maximum  thickness  of  8.84%  is  mounted  in  the  test  section  on  the  centerline 
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of  the  tunnel.  It  has  a  geometrical  angle  of  attack  of  1.18*’.  Upstream  there  are  two 
flappers  operating  at  the  reduced  frequency  k  =  3.6  and  with  the  amplitude  of  6  degrees. 
The  Reynolds  numbers  based  on  testfoil  chord  length  and  freestream  velocity  are  1.36  x  10® 
and  3.78  x  10®.  Free  stream  turbulence  intensity  was  measured  in  the  empty  test  section 
and  found  to  be  1%.  Transition  was  artificially  triggered  dX  xjC  =  0.105  on  both  suction 
and  pressure  sides  using  0.05  inch  diameter  epoxy  disks  of  0.008  inch  high  and  separated 
by  0.05  inch. 

Velocity  and  surface  pressure  were  measured  using  two-component  Laser  Doppler  Ve- 
locimeter  (LDV)  and  miniature  pressure  transducers.  For  the  unsteady  measurements, 
the  value  at  each  time  step  was  recorded  as  tan  average  over  250  periods.  There  are  180 
sample  points  over  one  period.  The  velocity  and  pressure  data  were  accumulated  over  a 
“bounding-box”  around  the  test-foil,  as  well  as  on  the  surface  and  at  certain  wake  locations. 
Friction  coefficients  Cj  data  was  deduced  from  surface  boundary  layer  measurements.  The 
standard  deviations  of  the  velocity  data  on  the  “bounding  box”  are  0.029  in  the  streamwise 
direction  and  0.016  in  the  vertical  direction.  The  unsteady  pressure  measurements  had  a 
standard  deviation  of  0.02  based  upon  absolute  pressure.  There  were  uncertainties  in  the 
unsteady  pressure  measurements  —  it  was  found  that  the  harmonic  content  of  the  surface 
pressure  responds  to  the  flapper  or  driving  frequency,  rather  than  to  the  hydrodynamic 
parameter  of  reduced  frequency.  It  is  very  possible  that  the  flappers  excited  a  mode  of 
tunnel  vibration  which  contaminated  the  higher  harmonics  of  pressure  measurements. 

We  have  performed  numerical  computations  for  both  steady  and  unsteady  flows  with 
Reynolds  number(3.78  x  10®).  The  experimental  setup  described  above  was  fully  simulated 
in  our  numerical  model  by  incorporation  of  the  data  provided  by  the  MIT  experimental 
group. 
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4.2  Computational  Grid,  Initial  and  Boundary  Con¬ 
ditions 


The  origin  of  the  coordinate  system  is  at  the  leading  edge  of  the  test  foil,  with  the 
a;— coordinate  positive  in  the  downstream  direction.  The  computational  domain  covers 
the  tunnel  with  the  inlet  positioned  upstream  of  the  leading  edge  dX  xjC  =  —0.261,  which 
coincides  with  the  location  where  experimental  measurement  of  the  flow  data  is  available, 
and  the  outlet  being  downstream  at  xjC  —  2.5.  Here,  C  is  the  hydrofoil  cord  length.  A 
C-type  grid  built  around  the  test  foil  uses  the  total  of  301x166  grid  points.  In  order  to 
assure  the  proper  boundary  layer  resolution  without  wasting  computational  resources  in 
the  far  flow  region,  the  grid  spacing  has  been  adjusted  by  hyperbolic  tangent  functions: 


where: 


and 


As(/)  = 


<l(I) 

A+{l-A),(/) 


-  I)} 

tanh(|) 


A  - 

^/Ai7 

As  a  result,  the  grid  spacing  near  the  surface  of  the  test  foil  is  much  denser,  with  approxi¬ 
mately  35  grid  points  in  the  boundary  layer.  It  should  be  mentioned  that  as  a  result  of  the 
this  adjustment  the  grid  spacing  at  the  tunnel  wall  is  too  coarse  to  resolve  any  boundary 
layer  features.  However,  at  these  high  Reynolds  number(fti  10®),  the  tunnel  wall  boundary 
layer  is  so  thin  that  there  is  no  influence  upon  the  flow  around  test  foil  which  is  of  interest 
here. 
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The  spatial  discretization  is  second-order  accuracy.  For  unsteady  flow,  the  time-stepping 
scheme  is  fully  implicit  first-order  Euler,  and  hence  absolutely  stable.  Assuming  that  qij 
is  a  time  dependent  scalar  quantity  at  the  {i,j)  grid  point,  the  time-discretization  is  done 
as  follows: 


Standard  non-slip  boundary  conditions  are  applied  on  the  foil  surface  and  tunnel  walls. 
Downstream  of  the  test  foil,  we  adopted  the  outlet  boundary  conditions  which  imply  zero 
normal  fluxes  of  all  flow  variables.  The  upstream  velocities  are  incorporated  from  the 
experiment  data.  This  data,  however,  is  available  only  in  the  section  between  y  jC  =  —0.306 
and  y/C  =  0.360  while  the  tunnel  walls  are  at  yjC  =  —0.525  and  y/C  —  0.525.  In  the 
steady  case,  the  experimental  data  was  interpreted  via  cubic  spline  to  computational  grid 
points.  In  the  region  between  tunnel  wall  and  the  flapper  wake  where  the  experimental 
data  is  not  available,  the  uniform  velocity  profiles  were  assigned.  In  the  case  of  unsteady 
flow,  the  MIT  FFX  experiment  provided  180  velocity  profile  samples  inside  the  flapper 
wake  at  the  upstream  boundary.  We  performed  a  harmonic  analysis  of  the  experimental 
data: 


U{x,y,t) 


20 
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and  extracted  20  harmonics  as  the  upstream  boundary  conditions.  In  the  region  between 
the  flapper  wake  and  the  tunnel  walls,  we  have  taken  two  different  extrapolations;  1)  simply 
setting  the  zeroth  harmonics  of  the  velocities  to  their  freestream  value  and  setting  all  the 
multiple  harmonics  to  zero,  and  2)  setting  the  zeroth  harmonics  to  some  freestream  value 
and  linearly  extrapolating  the  amplitude  and  the  phase  of  the  first  harmonics  from  the 
edge  of  the  flapper  wake  data: 


U{x,youter,t)  =  t/oo  +  Ui{x,yedge)  Sin{ujt  +  6  i{x ,  y  edge))  ■ 


Another  alternative  is: 

U{x,youter.t)  =  C/oo  ±  6^1  (o:,  +  0^{x,yedge)) 


^ y outer  5 


f>  /  X  sinh(  (2/i  =p  j/oufer))  '(  X  \  a  ( 
Vr{x,yedge) - ^Qshj'M) - sin(u;f  +  ^1  (a:,  yerf^e)) 


with  all  the  higher  harmonics  set  to  zero.  Finally,  as  for  the  steady  case,  profiles  of  all 
velocity  harmonics  were  projected  onto  the  computation  grid  using  cubic  spline  fit. 


4.3  Computational  Results 
4.3.1  steady  Flow 

Figures  4.1  to  4.3  show  the  computed  static  pressure  and  velocity  fields  respectively.  It 
takes  about  10,000  iterations  to  achieve  convergence  when  started  from  a  uniform  initial 
field.  We  discuss  here  various  quantitative  features  of  the  steady  flow. 

Test  box  data  and  computational  inlet  profiles. 

The  available  experimental  data  includes  average  velocities  at  several  stations  on  the 
sides  of  three  rectangular  nested  boxes  surrounding  the  test  foil.  In  each  direction,  the 
average  space  separation  between  the  probes  is  about  0.2%  ~  0.3%.  This  data  at  the 
upstream  face  of  the  box  was  used  to  set  up  the  inlet  velocity  conditions,  as  shown  in 
Figure  4.4.  The  experimental  data  corresponding  to  the  top,  bottom  and  the  downstream 
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faces  of  the  middle  box  were  used  to  compare  the  computed  mean  velocity  field  in  the  far 
field  region.  The  data  in  Figure  4.5  to  Figure  4.7  respectively  is  indicative  of  adequate  grid 
resolution  and  satisfactory  boundary  condition  setting. 

Wake  velocity  profiles. 

The  wake  velocity  profile  and  flow  inclination  plotted  in  Figure  4.8  demonstrate  good 
overall  agreement  with  the  experimental  data.  The  relative  flow  angle  computed  using 
RNG  K  —  e  model  clearly  shows  that  the  flow  direction  above  the  test  foil  is  not  aligned 
with  the  geometric  foil  surface.  In  particular,  the  flow  has  separated  in  a  certain  region  just 
upstream  of  the  trailing  edge.  This  pattern  is  confirmed  by  contour  plots  of  the  streamwise 
velocity  around  the  foil  trailing  edge  shown  in  Figure  4.11.  It  should  be  mentioned  that  the 
standard  K  —  e  model  completely  fails  to  predict  separation  near  the  trailing  edge,  which  is 
not  surprising  since  the  standard  K  —  t  model  tends  to  overpredict  eddy  viscosity  so  that 
the  turbulence  level  inside  the  boundary  layer  is  too  high  for  the  flow  to  separate.  We  plot 
in  Figure  4.12  the  turbulence  kinetic  energy  predicted  by  both  the  standard  K  —  t  model 
and  the  RNG  K  —  t  model.  Obviously,  the  standard  K  -  t  model  generates  much  stronger 
turbulence  energy  on  the  foil  surface  and  into  the  near  trailing  edge  wake  region. 

Surface  pressure. 

Figure  4.9  shows  the  surface  pressure  coefficients.  The  observed  variations  of  the  com¬ 
puted  Cp  curves  when  the  angle  of  attack  is  adjusted  between  1.18  and  1.34  degrees  is 
probably  due  to  the  uncertainty  in  the  experimental  conditions. 

Skin  friction  coefficients. 

Obviously,  prediction  of  the  skin  friction  coefficients  involves  adequate  quantitative  de¬ 
scription  of  the  entire  boundary  layer  and  transitional  region  thus  representing  the  ultimate 
test  of  a  transport  model  performance.  The  RNG  K  —  e  model  has  demonstrated  undis- 
putable  superiority  over  the  existing  transport  models. 

In  Figure  4.10,  we  plotted  the  comparison  of  skin  friction  coefficients  computed  using 
the  standard  and  the  RNG  K  —  e  model  against  experimental  measurements.  On  the 
suction  side  for  xjC  >  0.389,  where  experiment  data  were  available,  computation  with 
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the  RNG  K  —  e  model  produced  excellent  agreement  while  the  computation  results  using 
standard  K  —  e  model  severely  overpredicted  the  skin  friction.  The  results  from  the  RNG 
K  -  e  model  predicted  the  separation  point  (shown  by  the  location  where  skin  friction 
first  turns)  at  x/C  =  0.968,  whereas  the  experiment  data  showed  the  separation  point 
between  x/C  =  0.972  and  x/C  -  0.990.  On  the  other  hand,  standard  K  -  t  model  did 
not  indicate  any  separation  at  all  because  it  overproduces  turbulence  energy,  as  we  have 
discussed  above.  The  RNG  K  -  e  model  computation  also  predicted  transition  regions  at 
x/C  0.09  on  the  suction  side  and  ai  x/C  ct  0.13  on  the  pressure  side.  Since  the  MIT 
experiment  employed  transition-tripping  devices  on  the  testfoil,  rigorous  comparison  with 
the  experiment  data  is  not  possible. 

4.3.2  Unsteady  Flow 

The  period  of  flapper  motion  is  T  =  y  =  0.06217(s).  We  conducted  time-step  sensitivity 
study  and  found  that  ^  is  an  adequate  choice  for  the  time  step.  The  unsteady 

calculation  usually  started  from  a  converged  steady  field  and  it  would  take  typically  about 
3  ~  4  periods  to  reach  a  periodic  state.  For  the  purpose  of  time  series  analysis,  32  samples 
per  period  were  taken  in  most  of  our  numerical  simulations.  The  total  evolution  time  for 
the  time-dependent  structure  was  usually  4  periods. 

Zeroth  Harmonics 

We  compare  the  zeroth  harmonics  of  Cp  of  unsteady  flow  and  the  Cp  of  steady  flow 
in  Figure  4.13.  Our  computation  agrees  well  with  the  experimental  data,  except  for  the 
pressure  side  of  the  foil  where  there  is  5%  ~  10%  discrepancy  in  the  region  ^  >  0.6.  For  the 
most  part  of  the  foil,  the  steady  Cp  is  in  good  agreement  with  the  mean  Cp  of  the  unsteady 
flow.  Figure  4.14  shows  the  skin  friction  coefficients  from  steady  flow  computation  and  the 
zeroth  harmonics  of  the  skin  friction  coefficients  obtained  in  unsteady  flow  calculation.  We 
observe  that  the  averaged  C/  in  the  unsteady  flow  agrees  fairly  well  with  its  counterparts 
in  the  steady  flow  almost  everywhere  except  for  the  transition  region  where  the  mean  Cf  of 
unsteady  flow  shows  ~  15%  increase  over  the  Cj  of  steady  flow;  the  transition  point  on  the 
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suction  side  is  shifted  forward  in  the  streamwise  direction  by  about  5%  while  the  transition 
point  on  the  pressure  side  of  the  foil  almost  does  not  change.  Although  the  flapper  frequency 
[fj  z=  16)  is  close  to  the  natural  frequency  of  the  flow  ^  ~  14,  the  the  unsteadyness 
effect  is  minimal  due  the  fact  that  the  amplitude  of  the  unsteady  vertical  gust  on  the  inlet 
is  small  =  0.022).  Our  main  conclusion  is  that  the  steady  calculation  faithfully 

reproduces  the  average  quantities  in  the  unsteady  flow.  Although  further  experimental  and 
numerical  study  may  improve  the  accuracy  of  setting  inlet  conditions  in  unsteady  flows, 
we  believe  that  for  practical  applications,  steady  state  analysis  based  on  the  RNG  VLES 
model  will  provide  accurate  quantitative  description. 

Higher  harmonics 

We  compare  the  computed  first  harmonics  of  the  surface  pressure  coefficient  with  the 
experimental  data  in  Figure  4.15.  The  numerical  data  reproduces  well  the  qualitative 
behavior.  The  observed  10%  ~  20%  difference  between  the  computational  and  experimental 
data  for  the  harmonics  amplitude  is  largely  due  to  the  sensitivity  of  the  pressure  harmonics 
to  the  time  dependent  inlet  conditions.  Indeed,  after  using  a  different  interpolation  of 
the  inlet  data,  as  well  as  the  potential  model  approximation  for  the  region  between  the 
flapper  wake  and  the  tunnel  wall  where  experimental  data  is  not  available,  we  have  found 
that  the  pressure  harmonics  amplitudes  differ  by  the  same  10%  ~  20%.  Figure  4.16  shows 
the  amplitude  and  the  phase  of  the  skin  friction  coefficients  on  the  test  foil  surface.  The 
amplitude  shows  a  peak  near  the  leading  edge  between  ^  =  0.08  and  ^  =  0.16;  similar 
behavior  of  the  higher  harmonics  amplitude  has  been  confirmed  by  the  RNG  K  —  e  model 
computation.  Incidentally,  the  zeroth  harmonics  shows  that  the  unsteady  flow  undergoes 
transition  in  the  same  region. (  see  Figure  4.14  ). 


42 


.■r.".  mi 

illiS 


-5.00  -3.57  -2.14  -0.71  0.71  2.14  3.57  5.00 


Figure  4.1:  Steady  flow  static  pressure  contours  around  the  foil. 
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Figure  4.2;  Steady  flow  streamwise  velocity  contours. 
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Figure  4.3:  Steady  flow  vertical  velocity  contours. 
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Figure  4.6;  Velocity  profiles  at  the  bottom  face  of  the  test  box,  q  =  - 
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Figure  4.7:  Velocity  profiles  at  the  downstream  face  of  the  test  box, 
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Figure  4.8:  Wake  profiles  near  the  trailing  edge,  ^  =  1.05. 
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Figure  4.11:  Streamwise  velocity  contours  at  the  trailing  edgeof  the  foil,  Separation  region 
on  the  suction  side  is  clearly  visible. 
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Figure  4.12:  Turbulence  kinetic  energy  contours  at  the  trailing  edge  of  the  foil.  Red  color 
for  standard  K  —  t  model  indicates  a  level  of  ~  2%  —  5%  when  normalized  by  free  stream 
kinetic  energy,  while  red  color  for  RNG  K  —  e  model  represents  Ci;  2%.  Blue  color  shows 
~  10“^%  for  both  models. 
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Figure  4.13:  Zeroth  Harmonics  of  the  surface  pressure  coefficients. 


Figure  4.14:  Comparison  of  the  skin  friction  coefficient  of  steady  flow  and  the  mean  Cj  of 
unsteady  flow. 
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Figure  4.15:  First  Harmonics  of  the  surface  pressure  coefficients 
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